Update 150Nd to use the 3-stage Bateman equation for 151Sm levels - #130

diff --git a/periodictable/activation.py b/periodictable/activation.py
index b03fe95..394659c 100644
--- a/periodictable/activation.py
+++ b/periodictable/activation.py
@@ -227,6 +227,7 @@ Code original developed for spreadsheet by Les Slaback of NIST.
 
 
 from math import exp, log, expm1
+from sys import float_info
 import os
 from collections.abc import Callable, Sequence
 from typing import cast
@@ -336,6 +337,14 @@ class Sample:
         """
         After determining the activation, compute the number of hours required to achieve
         a total activation level after decay.
+
+        To simplify the decay time estimate, the 151Sm decay from 150Nd activation is
+        treated as 151Pm -> 151Sm instead of 151Nd -> 151Pm -> 151Sm, with the
+        missing 151Nd activity added to 151Sm at t=0 to compensate. This is
+        different from the decay calculation, which uses the full 3-stage Bateman
+        equation for the 151Sm activity. As a result, the decay time estimate will
+        be slightly long for short exposures and slightly short for long exposures,
+        depending on whether 151Sm activity exceeds the target threshold.
         """
         if not self.rest_times or not self.activity:
             return 0
@@ -354,13 +363,15 @@ class Sample:
         # dominate at long times, but at short times they will not affect the
         # derivative. Choosing a time that satisfies the longest decay time seems
         # to work well enough.
-        t_guess = -1e100  # Allow the initial guess to be negative
+        MIN_TIME = -1e100  # Allow the initial guess to be negative
+        t_guess = MIN_TIME
         initial_activity = 0.  # Accumulate the intensity at t=0.
+        activity_Nd151 = 0.  # Linter doesn't know activity_Nd151 is defined before use
 
         for k, (a, Ia) in enumerate(self.activity.items()):
             intensity = Ia[0]
             La = LN2/a.Thalf_hrs
-            # print(f"{k}: {a.daughter} {intensity=} {La=}")
+            # print(f"  {k}: {a.daughter} {intensity=} {a.Thalf_hrs=} {La=}")
             # Estimate activity at t=0, with a check that it hasn't decayed to zero.
             # The check is necessary because exp(λt) will exceed the floating point
             # number range when I(t) = I(0)exp(-λt) = 0.
@@ -371,21 +382,38 @@ class Sample:
                     Ip, Lp, _ = data[-1] # parent activity and decay constant
                     #print(f"{intensity=} {Ip=} {Lp=} {Ia[0]=} {La=}")
                     intensity -= Ip*(expm1((La-Lp)*t)/(1 - Lp/La)) if Ip > 0. else 0.
+            # Hack: Subtract scaled 151Nd activity from 151Sm so that we don't need
+            # to fit with the three stage Bateman equation. This works because there
+            # is a million-fold difference in half-lives. If there is enough activity
+            # in 151Nd to push 151Sm over the threshold then the 151Sm intensity will
+            # determine decay time, otherwise 151Pm will dominate. Because 151Sm activity
+            # is front-loaded, the integrated intensity for short times will be a little
+            # higher than expected, giving a slightly long decay time estimate when
+            # 151Sm activity is below the threshold.
+            if a.daughter.startswith("Nd-151"):  # Nd-151t
+                activity_Nd151 = intensity
+            elif a.isotope == "Nd-150" and a.daughter.startswith("Sm-151"): # Sm-151
+                # print(f"151Nd correction to 151Sm = {activity_Nd151 * a.Thalf_parent / a.Thalf_hrs}")
+                intensity += activity_Nd151 * a.Thalf_parent / a.Thalf_hrs
             initial_activity += intensity
             data.append((intensity, La, a.reaction))
 
-            # Daughter intensity will be transformed to granddaughter intensity
-            # over time. Assume it is instantaneous at t=0 for the purpose of guessing
+            # Daughter intensity with delayed beta decay will be transformed to granddaughter
+            # intensity over time. Assume it is instantaneous at t=0 for the purpose of guessing
             # the decay time of the granddaughter. The root finder will correct for
             # the discrepency.
-            I_buildup = data[-1][0]*La/data[-1][1] if a.reaction == "b" else 0.
-            decay_estimate = -log(target/(intensity + I_buildup))/La if intensity > 0. else 0.
-            t_guess = max(t_guess, decay_estimate)
-        #print("corrected at t=0", [Ia for Ia, La, reaction in data])
-
-        # No need to fit decay time if initial intensity is below target. Even with two-stage
-        # decay chains, we are only looking at daughter products with longer half life than
-        # the parent, so activity of the parent during burnup will decrease more than the
+            # Note: data[k-1] != data[-1] because we have already appended to data in the loop
+            # Note: "b" is never the first reaction, so data[k-1] is always valid
+            Ib = data[k-1][0]*La/data[k-1][1] if a.reaction == "b" else 0.
+            t_guess_k = -log(target/(intensity + Ib))/La if intensity + Ib > 0. else MIN_TIME
+            # print(f"  {t_guess_k=} {Ib=} {intensity+Ib=} {La=}")
+            t_guess = max(t_guess, t_guess_k)
+        # print("corrected at t=0", [Ia for Ia, La, reaction in data])
+        # print(f"{t_guess=}")
+
+        # No need to fit decay time if initial intensity is below target. Even with delayed
+        # beta decay, we are only looking at daughter products with longer half life than
+        # the parent, so beta decay of the parent decreases its activity more than the
         # the activity of the daughter can increase.
         if initial_activity <= target:
             return 0.
@@ -393,7 +421,7 @@ class Sample:
 
         # Build f(t) = total activity at time T minus target activity and its
         # derivative df/dt. f(t) will be zero when activity is at target.
-        # 2026-05-27 PAK: include post exposure granddaughter buildup
+        # 2026-05-27 PAK: include post exposure delayed beta activation
         def f(t):
             total = 0.
             for k, (Ia, La, reaction) in enumerate(data):
@@ -403,7 +431,7 @@ class Sample:
                     # and half-life from the previous row in the table.
                     Ip, Lp, _ = data[k-1]
                     intensity += Ip*(exp(-Lp*t) - exp(-La*t))/(1 - Lp/La)
-                # print(f" f{k}: {activity}")
+                # print(f" f{k}: {intensity}")
                 total += intensity
             # print(f"f({t}) = {total} - {target}\n")
             return total - target
@@ -417,7 +445,7 @@ class Sample:
                     Ip, Lp, _ = data[k-1]
                     delta += Ip*(La*exp(-La*t) - Lp*exp(-Lp*t))/(1 - Lp/La)
                 total += delta
-                # print(f" df{k}: {d_activity}")
+                # print(f" df{k}: {delta}")
             # print(f"dfdt({t}) = {total}\n")
             return total
         t, ft = find_root(t_guess, f, dfdt, tol=tol)
@@ -427,7 +455,7 @@ class Sample:
             #return 0. #  Decay time failed to compute; silently fail
             #return 1e100*365*24 # Return 1e100 rather than raising an error
             msg = (
-                f"Failed to compute decay time correctly ({percent_error:.1g} error)."
+                f"Failed to compute decay time correctly ({percent_error:.1g}% error)."
                 f" Please report material and activation parameters.")
             raise RuntimeError(msg)
 
@@ -669,15 +697,24 @@ def activity(
     rest period.
 
     Activations are listed in *isotope.neutron_activation*. Most of the
-    activations (n,g n,p n,a n,2n) define a single step process, where a neutron
-    is absorbed yielding the daughter and some prompt radiation. The daughter
+    activations (n,g n,p n,a) define a single step process, where a neutron is
+    absorbed yielding the daughter and some prompt radiation. The daughter
     itself will decay during exposure, yielding a balance between production and
     emission. Any exposure beyound about eight halflives will not increase
-    activity for that product.
+    activity for that product. Activity for daughter products may undergo
+    further neutron capture, reducing the activity of the daughter product but
+    introducing a grand daughter with its own activity (n,2n).
+
+    Active delayed beta products (b mode reactions) accumulate during
+    irradiation, eventually reaching equilibrium. Burn up, where the activated
+    product undergoes further capture before beta decay is not accounted for.
+    The delayed beta product increases in activity until the directly activated
+    product has decayed away following the 2-stage Bateman equation.
 
-    Activity for daughter products may undergo further neutron capture, reducing
-    the activity of the daughter product but introducing a grand daughter with
-    its own activity.
+    150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu has two transient daughters.
+    During irradiation, the activity for 151Sm is computed directly from 151Nd,
+    skipping 151Pm. After irradiation the 151Pm activity is subtracted from
+    151Sm, and decay calculated using the 3-stage Bateman equation.
 
     The data tables for activation are only precise to about three significant
     figures. Any changes to the calculations below this threshold, e.g., due to
@@ -718,9 +755,11 @@ def activity(
     if not hasattr(isotope, 'neutron_activation'):
         return result
 
-    # Hack to support b-mode 2-stage decay chain
-    # Relies on b line following directly after activation line
-    last_activity = 0.
+    # Hack: delayed beta calculation needs activity and decay rate
+    # of the parent. Since the "b" mode reaction line follows directly
+    # after the parent line, save the last activity and the last decay rate
+    # at each step of the loop. The 151Sm decay needs the 151Nd activity.
+    last_activity = last_lam = activity_151Nd = 0.
     for ai in isotope.neutron_activation:
         # Ignore fast neutron interactions if not using fast ratio
         if ai.fast and env.fast_ratio == 0:
@@ -784,7 +823,6 @@ def activity(
             #             = root * (lam*expm1(x2) - parent_lam*expm1(x1)) / (parent_lam - lam)
             # Checked for each b-mode production that small halflife results are
             # unchanged to four digits and Eu[151] => Gd[152] no longer fails.
-            # TODO: 150Nd => 151Nd => 151Pm => 151Sm => 151Eu
             activity = root/(parent_lam - lam) * (
                 lam*expm1(-parent_lam*exposure) - parent_lam*expm1(-lam*exposure))
             #print("N", parent_lam, "O", activity)
@@ -801,14 +839,17 @@ def activity(
             # Column S: resulting product [1/h]
             product_2n = lam if ai.reaction == '2n' else 0
             # Column T: activity if 2n mode
-            activity = root*lam*(parent_activity-parent_lam)*(
-                (exp(-lam_2n*exposure)
-                 / ((parent_activity-lam_2n)*(product_2n-lam_2n)))
-                + (exp(-parent_activity*exposure)
-                   / ((lam_2n-parent_activity)*(product_2n-parent_activity)))
-                + (exp(-product_2n*exposure)
-                   / ((lam_2n-product_2n)*(parent_activity-product_2n)))
-                )
+            #     activity = root*lam*(parent_activity-parent_lam)*(
+            #         (exp(-lam_2n*exposure)
+            #          / ((parent_activity-lam_2n)*(product_2n-lam_2n)))
+            #         + (exp(-parent_activity*exposure)
+            #            / ((lam_2n-parent_activity)*(product_2n-parent_activity)))
+            #         + (exp(-product_2n*exposure)
+            #            / ((lam_2n-product_2n)*(parent_activity-product_2n)))
+            #         )
+            # [CAOShurong 2026-10-15] Evaluate T as a three-rate convolution.
+            activity = root*lam*(parent_activity-parent_lam)*_three_rate_integral(
+                exposure, lam_2n, parent_activity, product_2n)
             #print("N", parent_lam, "P", effectiveXS, "Q", lam_2n,
             #      "R", parent_activity, "S", product_2n, "T", activity)
         else:
@@ -823,15 +864,24 @@ def activity(
             # Note: This cross-section always uses the total thermal flux
             effectiveXS = ai.thermalXS_parent + env.epithermal_reduction_factor*ai.resonance_parent
             # Column U: nv1s1t [neutrons]
-            U = (flux*3600)*(initialXS*1e-24)*exposure
+            #     U = (flux*3600)*(initialXS*1e-24)*exposure
             # Column V: nv2s2t+L2*t [neutrons]
-            V = ((env.fluence*3600)*(effectiveXS*1e-24)+lam)*exposure
+            #     V = ((env.fluence*3600)*(effectiveXS*1e-24)+lam)*exposure
             # Column W: L/(L-nvs1+nvs2) [unitless]
-            W = lam/(lam-(flux*3600)*(initialXS*1e-24)+(env.fluence*3600)*(effectiveXS*1e-24))
-            # Column X: W*(exp(-U)-exp(V)) if U,V > 1e-10 else W*(V-U+(V-U)*(V+U)/2)
-            # [PAK 2024-02-28] Rewrite the exponential difference using expm1()
-            #Xp = exp(-V)*expm1(V-U) if U>V else -exp(-U)*expm1(U-V)
-            X = W*exp(-V)*expm1(V-U) if U>V else -W*exp(-U)*expm1(U-V)
+            #     W = lam/(lam-(flux*3600)*(initialXS*1e-24)+(env.fluence*3600)*(effectiveXS*1e-24))
+            # Column X: W*(exp(-U)-exp(-V))
+            #     Original small-U,V form: W*(V-U+(V-U)*(V+U)/2)
+            # [PAK 2024-02-28] Rewrite the exponential difference using expm1():
+            #     X = W*exp(-V)*expm1(V-U) if U>V else -W*exp(-U)*expm1(U-V)
+            # [CAOShurong 2026-10-15] Rewrite U, V, W, X for higher accuracy.
+            # With a=production_rate, b=removal_rate, t=exposure:
+            #     U = a*t; V = b*t; W = lam/(b-a)
+            #     X = lam*(exp(-a*t)-exp(-b*t))/(b-a)
+            # The convolution also provides the finite equal-rate limit.
+            production_rate = (flux*3600)*(initialXS*1e-24)
+            removal_rate = (env.fluence*3600)*(effectiveXS*1e-24) + lam
+            # The exponential divided difference has a finite equal-rate limit.
+            X = lam*_two_rate_integral(exposure, production_rate, removal_rate)
             # Column Y: O if "b" else T if "2n" else L*X
             activity = root*X
             #print(f"{ai.isotope}=>{ai.daughter} {U=} {V=} {W=} {Xp=} {X=} {activity=} {lam=} {flux*initialXS*3600*1e-24+env.fluence*effectiveXS*3600*1e-24}")
@@ -847,30 +897,111 @@ def activity(
             #data = env.fluence, initialXS, flux, root, U, V, W, precision_correction
             #print(" ".join("%.5e"%v for v in data))
 
-        # TODO: 151Eu -> 152m2Eu -> 152Eu is missing b-mode 152Gd
-        # TODO: 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu is treated as a two stage decay
-        # It is present for 151Eu -> 152m1Eu isomer and 151Eu -> 152Eu ground state.
-        if ai.reaction == 'b':
+        # 2026-05-27 PAK: delayed beta decay such as 209Bi -> 210Bi -> 210Po -> 206Pb
+        # TODO: check 151Eu -> 152m2Eu -> 152Eu is missing b-mode 152Gd
+        # - It is present for 151Eu -> 152m1Eu isomer and 151Eu -> 152Eu ground state.
+
+        # print(f"{ai.isotope=} {ai.daughter=} {ai.reaction=}")
+        if ai.daughter.startswith("Nd-151"): # Nd-151t
+            # Track 151Nd activity at t=0 so we can correct 151Sm for short exposure
+            activity_151Nd = activity
+
+        if ai.isotope == "Nd-150" and ai.daughter.startswith("Sm-151"): # Sm-151
+            # Decay chain: 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu
+            # Because 151Sm uses 151Nd as its parent halflife, all the intensity from the
+            # transient 151Pm is accounted for in the 151Sm activity. We subtract the scaled 151Pm
+            # activity at t=0 from 151Sm. We use the three stage bateman equation to dribble activity
+            # from 151Nd and 151Pm into 151Sm.
+            # print(f"151Sm @ t=0 = {activity}")
+            # print(f"  - (activity_151Pm={last_activity})*{lam}/{last_lam} = {last_activity * lam/last_lam}")
+            # print(f"  = {activity - last_activity * lam/last_lam}")
+            activity -= last_activity * lam/last_lam
+
+            # Show the approximate activity in 151Sm after 151Nd and 151Pm have decayed
+            # print(f"151Sm @ t=280 h = {activity}")
+            # print(f"  + ({activity_151Nd=})*{lam}/{parent_lam} = {activity_151Nd * lam/parent_lam}")
+            # print(f"  + (activity_151Pm={last_activity})*{lam}/{last_lam} = {last_activity * lam/last_lam}")
+            # print(f"  = {activity + activity_151Nd * lam/parent_lam + last_activity * lam/last_lam}")
+
+            result[ai] = [
+                bateman3(Ti, activity, lam, last_activity, last_lam, activity_151Nd, parent_lam)
+                for Ti in rest_times
+            ]
+
+        elif ai.reaction == 'b':
             # Accumulate build up following the Bateman equation:
-            #   A_d(t) = λ_d/(λ_d - λ_p) A_p(0) (exp(-λ_p t) - exp(-λ_d t))
-            # Add this to decay of the granddaughter at the end of the exposure:
-            #          + A_d(0) exp(-λ_d t)
             result[ai] = [
-                activity*exp(-lam*Ti) + last_activity * (exp(-parent_lam*Ti) - exp(-lam*Ti)) / (1 - parent_lam/lam)
+                bateman2(Ti, activity, lam, last_activity, last_lam)
                 for Ti in rest_times
             ]
             # print(f"{ai.daughter} {activity=} {last_activity=}")
         else:
             result[ai] = [activity*exp(-lam*Ti) for Ti in rest_times]
-        # 2025-05-27 PAK: Hack to use 151Nd activation intensity rather than
-        # the 151Pm intensity when computing the buildup of the 151Sm granddaugter.
-        if ai.daughter != "Pm-151":
-            last_activity = activity
+        last_lam = lam  # For delayed beta computations (with 151Sm, last_lam != parent_lam)
+        last_activity = activity  # For delayed beta computations
         #print([(Ti, Ai) for Ti, Ai in zip(rest_times, result[ai])])
 
     return result
 
+def bateman2(t, act, lam, act_p, lam_p):
+    """
+    Accumulate build up following the Bateman equation::
+
+        A1(t) = A(0) exp(-λ t)
+        A2(t) = A_p(0) λ (exp(-λ_p t) - exp(-λ t))/(λ - λ_p)
+
+        A(t) = A1(t) + A2(t)
+
+    """
+    # simple decay
+    term1 = act*exp(-lam*t)
+    # two-stage decay
+    term2 = act_p * lam * _two_rate_integral(t, lam, lam_p)
+
+    return term1 + term2
+
+def bateman3(t, act, lam, act_p, lam_p, act_gp, lam_gp):
+    """
+    Accumulate build up following the three stage Bateman equation.
+    """
+    # simple decay
+    term1 = act * exp(-lam*t)
+    # two-stage decay
+    term2 = act_p * lam * _two_rate_integral(t, lam, lam_p)
+    # three-stage decay
+    term3 = act_gp * lam * lam_p * _three_rate_integral(t, lam, lam_p, lam_gp)
+
+    return term1 + term2 + term3
+
+def _two_rate_integral(t: float, a: float, b: float) -> float:
+    """Convolve exp(-a*t) and exp(-b*t), including coincident rates."""
+    if a > b:
+        a, b = b, a
+    x = (b-a)*t
+    # A zero or subnormal product loses precision before division by b-a.
+    if abs(x) < float_info.min:
+        return t*exp(-a*t)
+    # Divide first so a subnormal decay is not lost in the intermediate product.
+    return (-expm1(-x)/(b-a))*exp(-a*t)
+
+def _three_rate_integral(t: float, a: float, b: float, c: float) -> float:
+    """Convolve three exponentials without singular repeated-rate factors."""
+    a, b, c = sorted((a, b, c))
+    x, y = (b-a)*t, (c-a)*t
+    if y > 0.01:
+        # Use the largest separation in the denominator.
+        return (_two_rate_integral(t, a, b)-_two_rate_integral(t, b, c))/(c-a)
+
+    # Near the all-equal limit, expand the second divided difference.
+    # h_k(x,y) = sum(x**j*y**(k-j)); coefficient truncation error is < 3e-22.
+    coefficient, homogeneous, x_power, factorial = 0.5, 1., 1., 2
+    for order in range(1, 8):
+        x_power *= x
+        homogeneous = y*homogeneous + x_power
+        factorial *= order+2
+        coefficient += (-1)**order * homogeneous/factorial
+    return t*(t*exp(-a*t))*coefficient
+
 class ActivationResult:
     r"""
     *isotope* :